Liquid argon benchmarks

Sebastian Micluța-Câmpeanu, Mikhail Vaganov

The purpose of these benchmarks is to compare several integrators for use in molecular dynamics simulation. We will use a simulation of liquid argon form the examples of NBodySimulator as test case.

using ProgressLogging
using NBodySimulator, OrdinaryDiffEq, StaticArrays
using Plots, DataFrames, StatsPlots

function setup(t)
    T = 120.0 # K
    kb = 1.38e-23 # J/K
    ϵ = T * kb # J
    σ = 3.4e-10 # m
    ρ = 1374 # kg/m^3
    m = 39.95 * 1.6747 * 1e-27 # kg
    N = 216
    L = (m*N/ρ)^(1/3)
    R = 2.25σ
    v_dev = sqrt(kb * T / m) # m/s

    _L = L / σ
     = 1.0
     = 1.0
    _m = 1.0
    _v = v_dev / sqrt(ϵ / m)
    _R = R / σ

    bodies = generate_bodies_in_cell_nodes(N, _m, _v, _L)
    lj_parameters = LennardJonesParameters(, , _R)
    pbc = CubicPeriodicBoundaryConditions(_L)
    lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => lj_parameters));
    simulation = NBodySimulation(lj_system, (0.0, t), pbc, /T)

    return simulation
end
setup (generic function with 1 method)

In order to compare different integrating methods we will consider a fixed simulation time and change the timestep (or tolerances in the case of adaptive methods).

function benchmark(energyerr, rts, bytes, allocs, nt, nf, t, configs)
    simulation = setup(t)
    prob = SecondOrderODEProblem(simulation)
    for config in configs
        alg = config.alg
        sol, rt, b, gc, memalloc = @timed solve(prob, alg();
            save_everystep=false, progress=true, progress_name="$alg", config...)
        result = NBodySimulator.SimulationResult(sol, simulation)
        ΔE = total_energy(result, t) - total_energy(result, 0)
        energyerr[alg] = ΔE
        rts[alg] = rt
        bytes[alg] = b
        allocs[alg] = memalloc
        nt[alg] = sol.destats.naccept
        nf[alg] = sol.destats.nf + sol.destats.nf2
    end
end

function run_benchmark!(results, t, integrators, tol...; c=ones(length(integrators)))
    @progress "Benchmark at t=$t" for τ in zip(tol...)
        runtime = Dict()
        ΔE = Dict()
        nt = Dict()
        nf = Dict()
        b = Dict()
        allocs = Dict()
        cfg = config(integrators, c, τ...)

        GC.gc()
        benchmark(ΔE, runtime, b, allocs, nt, nf, t, cfg)
        get_tol(idx) = haskey(cfg[idx], :dt) ? cfg[idx].dt : (cfg[idx].abstol, cfg[idx].rtol)

        for (idx,i) in enumerate(integrators)
            push!(results, [string(i), runtime[i], get_tol(idx)..., abs(ΔE[i]), nt[i], nf[i], c[idx]])
        end
    end
    return results
end
run_benchmark! (generic function with 1 method)

We will consider symplectic integrators first

symplectic_integrators = [
    VelocityVerlet,
    # VerletLeapfrog,
    PseudoVerletLeapfrog,
    McAte2,
    # CalvoSanz4,
    # McAte5,
    Yoshida6,
    KahanLi8,
    # SofSpa10
]
5-element Array{DataType,1}:
 VelocityVerlet
 PseudoVerletLeapfrog
 McAte2
 Yoshida6
 KahanLi8

Let us run a short simulation to see the cost per timestep for each method

config(integrators, c, τ) = [ (alg=a, dt=τ*cₐ) for (a,cₐ) in zip(integrators, c)]

t = 35.0
τs = 1e-3

# warmup
c_symplectic = ones(length(symplectic_integrators))
benchmark(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), 10.,
    config(symplectic_integrators, c_symplectic, τs))

results = DataFrame(:integrator=>String[], :runtime=>Float64[], =>Float64[],
    :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
run_benchmark!(results, t, symplectic_integrators, τs)

5 rows × 7 columns

integratorruntimeτEnergyErrortimestepsf_evalscost
StringFloat64Float64Float64Int64Int64Float64
1VelocityVerlet58.40980.0010.0036656235000700021.0
2PseudoVerletLeapfrog59.21550.0010.0208884350001050021.0
3McAte257.67130.0010.00980374350001050021.0
4Yoshida6229.8960.0010.0219439350005250021.0
5KahanLi8530.5120.0010.02239233500012250021.0

The cost of a timestep can be computed as follows

c_symplectic .= results[!, :runtime] ./ results[!, :timesteps]
c_Verlet = c_symplectic[1]
c_symplectic /= c_Verlet
5-element Array{Float64,1}:
 1.0
 1.0137927860016844
 0.9873557593070432
 3.935903166715239
 9.082574686069826

were we have normalized the cost to the cost of a VelocityVerlet step.

Let us now benchmark the solvers for a fixed simulation time and variable timestep

t = 40.0
τs = 10 .^range(-4, -3, length=5)

results = DataFrame(:integrator=>String[], :runtime=>Float64[], =>Float64[],
    :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
run_benchmark!(results, t, symplectic_integrators, τs, c=c_symplectic)

25 rows × 7 columns

integratorruntimeτEnergyErrortimestepsf_evalscost
StringFloat64Float64Float64Int64Int64Float64
1VelocityVerlet679.9420.00010.002797314000008000021.0
2PseudoVerletLeapfrog671.5620.0001013790.00030127439455811836761.01379
3McAte2683.1449.87356e-50.00050904540512312153710.987356
4Yoshida6684.5650.000393590.012887410162915244373.9359
5KahanLi8665.0120.0009082570.03578514404115414379.08257
6VelocityVerlet378.8160.0001778280.003945722249374498761.0
7PseudoVerletLeapfrog375.1860.0001802810.001470332218776656331.01379
8McAte2386.2710.0001755790.001859512278186834560.987356
9Yoshida6382.2060.0006999140.0289537571508572523.9359
10KahanLi8381.6620.001615140.00651373247668668129.08257
11VelocityVerlet212.6690.0003162280.01342051264922529861.0
12PseudoVerletLeapfrog214.3160.0003205890.01133411247713743151.01379
13McAte2217.3650.0003122290.003279971281113843350.987356
14Yoshida6218.7570.001244640.013221321384820723.9359
15KahanLi8209.7680.002872160.0680963139274874479.08257
16VelocityVerlet121.4050.0005623410.00608374711321422661.0
17PseudoVerletLeapfrog119.330.0005700980.00520789701642104941.01379
18McAte2124.6520.0005552310.00932556720432161310.987356
19Yoshida6121.9410.002213320.000567354180732710973.9359
20KahanLi8119.6320.005107510.18143178322741229.08257
21VelocityVerlet67.85050.0010.011078340001800041.0
22PseudoVerletLeapfrog66.92390.001013790.0224749394561183701.01379
23McAte268.74520.0009873560.00814607405131215410.987356
24Yoshida669.4610.00393590.0120663101631524473.9359
25KahanLi867.66910.009082570.30801644051541779.08257

The energy error as a function of runtime is given by

@df results plot(:EnergyError, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")

Looking at the runtime as a function of timesteps, we can observe that we have a linear dependency for each method, and the slope is the previously computed cost per step.

@df results plot(:timesteps, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Number of timesteps", ylabel="Runtime (s)")

We can also consider a longer simulation time

t = 100.0

τs = 10 .^range(-4, -3, length=5)

results = DataFrame(:integrator=>String[], :runtime=>Float64[], =>Float64[],
    :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
#run_benchmark!(results, t, symplectic_integrators, τs, c=c_symplectic)

The energy error as a function of runtime is given by

#@df results plot(:EnergyError, :runtime, group=:integrator,
#    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")

We can also look at the energy error history

function benchmark(energyerr, rts, ts, t, configs)
    simulation = setup(t)
    prob = SecondOrderODEProblem(simulation)
    for config in configs
        alg = config.alg
        sol, rt = @timed solve(prob, alg(); progress=true, progress_name="$alg", config...)
        result = NBodySimulator.SimulationResult(sol, simulation)
        ΔE(t) = total_energy(result, t) - total_energy(result, 0)
        energyerr[alg] = [ΔE(t) for t in sol.t[2:end]]
        rts[alg] = rt
        ts[alg] = sol.t[2:end]
    end
end

ΔE = Dict()
rt = Dict()
ts = Dict()
configs = config(symplectic_integrators, c_symplectic, 2.3e-4)
benchmark(ΔE, rt, ts, 45., configs)

plt = plot(xlabel="Rescaled Time", ylabel="Energy error", legend=:topleft);
for c in configs
    plot!(plt, ts[c.alg], abs.(ΔE[c.alg]), label="$(c.alg), $(rt[c.alg])s", yscale=:log10)
end
plt

Now, let us compare some adaptive methods

adaptive_integrators=[
    # DPRKN
    DPRKN6,
    # DPRKN8,
    # DPRKN12,
    # others
    Tsit5,
    # Vern7,
    # Vern9
]

config(integrators, c, at, rt) = [ (alg=a, abstol=at, rtol=rt) for a in integrators]

t = 35.0
ats = 10 .^range(-20, -14, length=10)
rts = 10 .^range(-20, -14, length=10)

# warmup
benchmark(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), 10.,
    config(adaptive_integrators, 1, ats[1], rts[1]))

results = DataFrame(:integrator=>String[], :runtime=>Float64[], :abstol=>Float64[],
    :reltol=>Float64[], :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
run_benchmark!(results, t, adaptive_integrators, ats, rts)

20 rows × 8 columns

integratorruntimeabstolreltolEnergyErrortimestepsf_evalscost
StringFloat64Float64Float64Float64Int64Int64Float64
1DPRKN632.90541.0e-201.0e-200.03639266035452381.0
2Tsit518.19851.0e-201.0e-200.14443281217171.0
3DPRKN632.9954.64159e-204.64159e-200.002142196064456931.0
4Tsit518.43534.64159e-204.64159e-200.2873783262216871.0
5DPRKN631.65422.15443e-192.15443e-190.06185275916442371.0
6Tsit518.84752.15443e-192.15443e-190.9426053317222811.0
7DPRKN631.99171.0e-181.0e-180.0122435948445731.0
8Tsit517.9591.0e-181.0e-180.7980383231214771.0
9DPRKN632.14734.64159e-184.64159e-180.03026645992449511.0
10Tsit518.44084.64159e-184.64159e-180.9951483236215311.0
11DPRKN631.51062.15443e-172.15443e-170.0567565894440621.0
12Tsit518.5042.15443e-172.15443e-170.7323073281219571.0
13DPRKN630.49191.0e-161.0e-160.05827335767426901.0
14Tsit518.08881.0e-161.0e-160.8056523270216871.0
15DPRKN632.92334.64159e-164.64159e-160.03697945928444751.0
16Tsit518.22094.64159e-164.64159e-160.2463023256216211.0
17DPRKN631.9282.15443e-152.15443e-150.02091995917440761.0
18Tsit518.63432.15443e-152.15443e-150.02145933269217291.0
19DPRKN631.40971.0e-141.0e-140.1000835872437681.0
20Tsit518.37411.0e-141.0e-140.3145473259217831.0

The energy error as a function of runtime is given by

@df results plot(:EnergyError, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")

If we consider the number of function evaluations instead, we obtain

@df results plot(:EnergyError, :f_evals, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Number of f evals")

We can also consider a longer simulation time

t = 100.0

ats = 10 .^range(-20, -14, length=10)
rts = 10 .^range(-20, -14, length=10)

results = DataFrame(:integrator=>String[], :runtime=>Float64[], :abstol=>Float64[],
    :reltol=>Float64[], :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
#run_benchmark!(results, t, integrators, ats, rts)

The energy error as a function of runtime is given by

#@df results plot(:EnergyError, :runtime, group=:integrator,
#    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")